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Abstract 

Non-genetic perturbations, such as environmental change or developmental noise, can induce novel phenotypes. 
If an induced phenotype confers a fitness advantage, selection may promote its genetic stabilization. Non-genetic 
perturbations can thus initiate evolutionary innovation. Genetic variation that is not usually phenotypically visible 
may play an important role in this process. Populations under stabilizing selection on a phenotype that is robust to 
mutations can accumulate such variation. After non-genetic perturbations, this variation can become a source of 
new phenotypes. We here study the relationship between a phenotype's robustness to mutations and a population's 
potential to generate novel phenotypic variation. To this end, we use a well-studied model of transcriptional regu- 
lation circuits. Such circuits are important in many evolutionary innovations. We find that phenotypic robustness 
promotes phenotypic variability in response to non-genetic perturbations, but not in response to mutation. Our 
work suggests that non-genetic perturbations may initiate innovation more frequently in mutationally robust gene 
expression traits. 



Introduction 

Two main perspectives exist about the origin of evolu- 
tionary innovations. The orthodox "genotype-first" per- 
spective emphasizes the role of mutations in the pro- 
duction of new phenotypes. In this perspective, muta- 
tions produce individuals with novel phenotypes whose 
frequency in a population may increase through natu- 
ral s election. Th e heterodox "phenotype-first" p e rspec- 
tive dHal l200il iMoczekl 120071 iNewman et all | 2006[ 
Palmed 12004c IPigliucci et all 120061 IPrice et all 12003c 
West-Eberhardi [l989l 120031) emphasizes the role of non- 
genetic perturbations, such as exposure to different tem- 
peratures, diets, or biotic interactions. Non-genetic 
perturbations also comprise fluctuations in an organ- 
ism's internal "microenvironment", such a s gene activity 
changes caused by noisy gene expression (Elowit zet all 
2002; Mc Adams and Arkinlll997l:lRai et alll20ld) . 



The phenotype-first perspective is based on the ob- 
servation that organisms often have highly plastic phe- 
notypes. That is, the same genotype has the potential to 
produce different phenotypes depending on non-genetic 
influences. Thus, a perturbation can trigger a plastic phe- 
notypic response in some individuals of a population. If 
the resulting novel phenotype provides a benefit to its 
carrier, it facilitates survival. Then, selection may in- 
crease the frequency of new or already existing genetic 



variants that exaggerate, refine or "stabilize" this phe- 
notype by making it independent of non-genetic factors. 
Waddington coined t he term genetic assi milation for this 
stabilization process ( 1 Waddin gtonl 1 1 95 3b . 



Increasing amounts of evidence support the impor 
tance for innovation of traits induce d by non-genetic fac 
tors ( f or di s sentin g opinions, see 
12003 



53, Il999h . 



de Jong and Croziei 



that assimilation can 



First, theoretical work shows 
occur under broad condi- 



tions ( Ciliberti et al. 


2007bt Espinosa-Soto et al., 


201C ; Lande. 2009 


; Masel, 20041 Ricel 1998: 


Siegal and Bergman. 2002; Wagnen, 1996c Wagner et al., 


1997). Second, laboratory evolution experiments 



show that assimilation d oes occur (Eldaretal., 2009; 



Rutherford and Lindquist, Il998t ISuzuki and Niihout , 
2006c IWaddingtonl Il953l Il956l) . Third, studies in natu- 
ral populations suggest that genetic assim ilation of traits 
induced by no n- genetic factors is not rare (Palmer, 2004; 
Pigliucci and Murrenl 120031 IWest-Eberhardl \2003 . For 
example, taxa with genetically determined dextral or 
sinistral morphologies are frequently derived from taxa 
in which the direction of the asymmetry is no t geneti- 
cally fixed, but where it is a plastic response dPalmei , 
19961 120041) . This occurs for many traits, such as the 
side on which the eye occurs in flat fishes (Pleuronecti- 
formes), and the side of the larger first claw in decapods 
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(Thalassinidea) (Palmer, 119961) . Transitions like these 
indicate genetic assimilation of a direction of asymme- 
try originally induced by non-inheritable factors. More 
generally, traits where fixed differences among closely 
related species are mirrored by plastic variation within 
populations are good candidates for genetic assimilation. 
For example, amph i bian tr aits, such as gut morphology 
( Ledon-Rettig et all 2008 ) and limb lengt h and snout 
length (Go mez-Mestre and Buchholz , 2006), follow this 
pattern. 

A system is robust to genetic or non-genetic 
perturbations if its phenotype does not change when 
perturbed. Mutational robustness and robustness 
to non-genetic perturbation s are correlated w it h one 
another in m any cases ( Ancel and Fon tanal 2000; 
Ciliberti et all l2007bt Ide Visser et all 120031: iLehnerl 
20 id iMeikleiohn and HartlL 120021: IProulx et all l2007t 
Rutherford and Lindauist|. 1 19981). although exc e ptions 



exist ^Cooper et all 120061 : 
MaselandSiegalil2009b . 



Fraser and Schadt , 2010t 



The ability to produce evolutionary innovation, 
linked to the robustness of a biological sys - 



tem dAncel and Fontanal l2000t Iciliberti et all 12007a : 



Draghi and WagnerL 120091: IWagnerL 120051 l2008al) . At 
first sight, robustness seems to hamper innovation. First, 
a mutationally robust system produces less phenotypic 
variation in response to mutations . It may thus not facil- 
itate the genotype-first sce nario (ICiliberti et all l2007at 
Draghi and WagnerL 120091) . Second, a system robust 
to non-genetic factors shows little phenotypic plasticity. 
Thus, it may not support innovation under the phenotype- 
first scenario. However, the role of robustness in innova- 
tion is subtler than it may seem. This becomes evident 
when one considers how genotypes and their phenotypes 
are organized in a space of genotypes. 

Genotypes exist in a vast space of possible geno- 
types. Two genotypes are neighbors in this space if 
one can be transformed into the other by a single mu- 
tation. The distribution of phenotypes in genotype space 
shows some qualitative similarities for different kinds of 
systems, from RNA and protein molecules to metabolic 
networks and transcriptional regulation circuits. First, 
in these systems large sets of genotypes produce the 
same phenotype. Each of these sets can be traversed 
through single mutation steps that leave the phenotype 
unchanged. Such a set is also referred to as a neu- 



tral network or genotype network dSchuster et aUll994l) . 



Second, different neighborhoods of the same genotype 
netwo rk contain genotypes with very different geno- 
types dCiliberti et all l2007at iFerrada and WagnerL 1 20081: 
Lipman and Wilbur, 1991; Rodrigues and Wagner, 2009; 
Schult es and Bartelll2000tlSchuster et allll994llWagnerl 
20083) - 

To understand how mutational robustness relates to 
a system's ability to produce evolutionary innovations, 
it is useful to distinguish between the mutational robust- 



ness of a genotype and that of a phenotype. A geno- 
type G\ is mutationally more robust than another geno- 
type G2, if G± is more likely to maintain the same phe- 
notype than G2 in response to mutation. By extension, 
a phenotype Pi is mutationally more robust than P2 if 
the genotypes that produce Pi preserve Pi, on average, 
more often than the genotypes adopting P2 preserve P2 
in response to mutations. Perhaps surprisingly, muta- 
tional phenotypic robustness can facilitate the produc- 
tion of novel ph enotypes for RNA structure phenotypes 
(Wagner, 2008b). The reason is that genotypes with a 
more robust phenotype form larger genotype networks 
and have, on average, more neighbors with the same 
phenotype. A population of such genotypes encoun- 
ters relatively few deleterious mutations that would slow 
its diversification and spreading through genotype space 
(while preserving its phenotype). It thus attains a higher 
genotypic diversity which translates into greater pheno- 
typic variability in response to mutations, even though 
every single genotype may have access to fewer other 



every single genotype may 
phenotypes dWagne r. 2008a) 



This mechanism, although c orroborated for RNA 
and protein structural phenotypes ( Ferrada and Wagner . 
20081 IWagnerl l2008bl) may not lead to increased phe- 



notypic variability in all systems. The reason is that it 
depends on how many different and unique phenotypes 
the neighborhood of different genotypes contains, and on 
how rapidly populations can spread through a genotype 
network. In other words, it depends on the organization 
of genotype networks in genotype space, which may dif- 
fer among different system classes. For example, a re- 
cent theoretical analysis suggests that the relationship be- 
tween phenotypic mutational robustness and the potential 
to generate phenotypic variation through m utation need 
not even be monotonic ( Drag hi et al. L Eoioh . 

The above considerations pertain to phenotypic 
variability in response to mutations. One might think 
that phenotypic variability in response to non-genetic 
perturbations may behave similarly since robustness to 
mutations a nd to non-genetic fac t ors are often positively 
correla ted ( Ancel and Fontanal 2000; Ciliberti et al. 



2007b; 



Ide Visser et alJ 



2003 



Lehnerl 



2010; 



Meikleiohn and HartlL 120021 iRutherford and Lm dquist. 
19981) . However, we show that this is not necessarily 



so for transcriptional regulation circuits, which exist on 
a higher level of organization than individual evolving 
molecules. Such circuits direct the production of specific 
gene activity patterns at particular times and places in 
the developing organism. Changes in the expression of 
their genes are involved in many evolutionary innova- 
tions dDavidson and Erwin , 20061: Shubin et al. , 2009). 



We study a generic computational model of transcrip- 
tional regulation in which the genotypes correspond to 
the cis-regulatory interactions in a transcriptional circuit. 
The phenotypes correspond to the gene activity pattern a 
circuit produces. 
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For this system, we show that high phenotypic ro- 
bustness to mutations increases the number of novel 
expression phenotypes a circuit can produce in re- 
sponse to non-genetic perturbations. Thus, pheno- 
typic robustness to mutation facilitates innovation un- 
der the phenotype-first scenario. It does so by al- 
lowing the accumulation of genetic variation that is 
not observed phenotypically under typical conditions, 
but th at may be exposed after non-genetic perturba- 
tions (de Visser et all 2003t iGibson and Dworkin , 2004 ; 
Masej l2006t iMasel and Siegali 120091) . 



Methods 
Model 

The model represents a regulatory circuit of N genes, 
where each gene's activity is regulated by other genes 
in the circuit. The circuit's genotype is defined by a 
real-valued matrix A = (ay), in which non-zero el- 
ements represent regulatory interactions between genes 
(Fig. la). An interaction (ay ^ 0) means that the ac- 
tivity of gene j can either have a positive (ay > 0) or a 
negative (ay < 0) effect on the activity of gene i. We use 
m to refer to the number of interactions in a given circuit, 
and c to its interaction density, i.e. to the number of inter- 
actions m divided by the maximum possible number of 
interactions TV 2 . A vector St — (sj , s\ N ') describes 
the activity state of the circuit at time t. 

The activity of the genes in the circuit changes ac- 
cording to the difference equation 



,(*) 

S t + T 



N 



(1) 



where a(x) equals -1 when x < 0, it equals 1 when 
x > 0, and it equals when x = 0. 

Despite its level of abstraction, variants of this 
model have proven useful for studying the evolution 
of robustness in gene regulatory c i rcuits dCiliberti et all 
2007bt iMartin and Wagnerl l2008t ISiegal and Bergmanl 



2002 



Wagner, 



1996), the effect of recombination on 



the p roduction of negat i ve ep istasis (lAzevedo et al 



2006; IMartin and Wag ner. 2009), t he evolution of mod- 



ularit y in gene circuits dEspinosa-Soto and Wagnei , 
20101) and the evolution of new gene activity pat - 
terns dCiliberti et all l2007al braghi and Wagnerl 120091) . 



Similar models have also been successfully used to 
predict the dynami cs of developmental processes in 
plants and anim als (Mend oza and Alvarez-Buvllaill998t 



Miolsness et allll991l) . 

We consider circuits that start their dynamics from 
a particular initial gene expression state so- One can 
view this initial state as being specified by factors exter- 
nal to the circuit, be they environmental factors, signals 
from adjacent cells, maternal regulators, or any genes 
"upstream" of the circuit. The phenotype is the stable 



gene activity pattern that a circuit attains when start- 
ing from so- We here focus on circuits producing such 
stable patterns and disregard circuits with more complex 
dynamics. We define a circuit's fitness by a function that 
increases steeply with the similarity of its phenotype Soo 
to a reference activity state s^*, that is considered the 
optimal phenotype in a given environment. We consider 
circuits that attain the same Soo as equal with respect to 
their gene expression phenotype. Under this assumption, 
a mutation that transforms two such circuits into one an- 
other would be neutral with respect to this phenotype 
(Fig. lb). 

Determination of 1 -mutant neighborhoods 

In several of our analyses, we explored properties of the 
circuits that differ from a reference circuit genotype G 
by one single mutation. For each entry ay in the matrix 
A of G we considered the following cases: i) if ay = 
we check the phenotype of two mutants, one in which 
ay < 0, and one in which ay > 0; ii) if ay ^ we 
also check the phenotype of two mutants, one in which 
an interaction is lost (ay = 0), and one in which we 
change the value of ay while keeping its sign unchanged. 
Among all the variants in the one-mutation neighbor- 
hood, we allowed exclusively those that maintained the 
number of interactions within an interval [m_ , m+], thus 
keeping interaction density at a value close to c. In all 
cases m + — m_ = 5. We discarded circuits that did 
not attain a steady-state gene activity pattern in this ap- 
proach. Whenever a new non-zero value was required for 
a given ay, we chose an N(0,1) pseudorandom number, 
and forced its sign if needed. We defined the robustness 
to mutations of a genotype G as the fraction of G's 1- 
mutant neighbors that produce the same phenotype as G 
when their dynamics start from the initial state sq- To 
assess the phenotypes that G can access through muta- 
tions we registered and counted all the different pheno- 
types produced by the set of single mutant circuits that 
neighbor the reference circuit G. The approach is read- 
ily extended to entire populations, where we determined 
all the different phenotypes that occur in the neighbor- 
hood circuits in the population. We counted phenotypes 
that occurred in the neighborhood of two or more circuits 
only once. 

Evolving populations 

For the model we use, a given pair of initial and final ex- 
pression states (sq, s^*) is representative of all pairs with 
the same fraction d of individual genes' expression val- 
ues that differ between sq and s^*. For a pre-specified 
d, we thus chose an arbitrary su ch pair, and followed p re- 
viously established procedures dCiliberti et al. , 2007bl) to 
identify a circuit genotype G that is able to drive the sys- 
tem from so to The regulatory interactions in the 
initial genotype G are real numbers sampled from a nor- 
mal distribution with mean and standard deviation 1, 
i.e., an N(0,1) distribution. After having identified one 
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such genotype G, we created a population of 200 copies 
of it, and subjected this population to repeated cycles 
("generations") of mutations (with a probability of mu- 
tation of fi = 0.5 per circuit), and strong stabilizing se- 
lection on s£§* . Whenever a circuit underwent mutation, 
we picked one of the circuit's 1 -mutation neighbors at 
random (see above). We allowed exclusively those mu- 
tations that maintained the number of interactions within 
a predetermined interval [m_,m+], thus keeping inter- 
action density around a predetermined fraction c of non- 
zero interactions among all possible (N 2 ) interactions. 
In this study, m + — m_ = 5 in all cases. 

Throughout, we interpret a circuit's "fitness" as a 
survival probability. We followed the regulatory dynam- 
ics of each gene circuit with s n as init ial condition. As in 
previous work dCiliberti et all 12007 a b). we disregarded 
genotypes that did not produce fixed-point equilibrium 
states, or that produced phenotypes in which the activity 
of a gene was equal to zero (neither active nor inactive). 
We assigned genotypes that attained the reference pattern 
s£g* a maximally possible fitness of 1 . Thus, s^* repre- 
sents a pre-determined optimal gene expression state. We 
assigned genotypes that attained an equilibrium state Sqq 
that differed from s£g* in the activity state of k genes a 
fitness equal to (1 — k/N) 5 , which ensures a steep de- 
crease in survival probability even for small deviations 
from s£g*. Each generation, we constructed a new pop- 
ulation by sampling individuals with replacement from 
the previous generation, and subjecting copies of them 
to mutation with a probability /i. We kept each of these 
new individuals with a probability equal to its fitness, and 
continued sampling until the newly generated population 
had 200 members. For all the populations we study, we 
let the initial population of identical genotypes evolve for 
10 4 generations under selection for s£g', before collect- 
ing any simulation data. This allows the population to 
erase any traces of the initial genotype and to reach a 
plateau where phenotypic variability in response to ei- 
ther mutations or non-genetic perturbations varies little 
across generations. 

We define the genotypic distance between two cir- 
cuits as the minimum number of mutations needed to 
transform one circuit into the other, normalized by the 
maximally possible number of such mutations for cir- 
cuits with the same number of interactions. 

Noisy dynamics 

For each circuit in a population, we generated 5N dy- 
namic trajectories, each of which started from so. For 
each of these trajectories, and for each step of the regu- 
latory dynamics, we perturbed the activity of a randomly 
picked gene with a probability of 0.5. We then followed 
each trajectory until an activity pattern s had consecu- 
tively repeated itself, and labeled this pattern as s^. We 
then counted the number of different fixed-point equi- 
librium states that each circuit could attain in these 5iV 



trajectories. 
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Fig. 1. Gene regulatory circuit model, (a) A gene regulatory circuit. 
Black bars indicate genes that encode proteins which regulate the 
activity of other genes in a hypothetical circuit. The regulatory 
interactions are described by a matrix A = (a^j ). An interaction 
means that the activity of gene j can either have a positive (ciij > 0, 
red rectangles) or a negative (a^j < 0, blue rectangles) effect on the 
activity of gene i. (b) Gene circuits that differ in a single interaction 
are neighbors in genotype space. . Each large circle surrounds a 
distinct gene regulatory circuit. Red arrows represent activating 
interactions, and blue lines represent repressing interactions between 
different genes (black rectangles). Dashed lines represent the 
interactions that are necessary to convert the indicated circuits into the 
middle circuit, a) and b) are modified with permission from 
(Ciliberti et al., 2007b). (c) Example of novel phenotypes caused by 
three kinds of perturbations. A reference gene circuit produces 
phenotype Soo, in which four genes are active (genes 5-8; yellow) and 
four genes are inactive (genes 1-4; red). In a one-mutation 
neighborhood (all circuits that differ from the reference one by a 
single interaction) we find four phenotypes sJJ? 1 " different from the 
original Sac (left). If we perturb the system state of the reference 
circuit without altering its genotype, other novel phenotypes are 
encountered (center and right panels). The perturbations we used are 
either all single-gene perturbations in the initial condition so (center), 
or perturbations of the dynamical trajectory of the circuit ('noisy 
dynamics'; right). 

Random sampling of genotypes in genotype networks 

In order to sample properties of a given genotype net- 
work uniformly, we performed a random mutational walk 
restricted to this genotype network, that is, to circuits that 
attain a given s^J* from the initial state sq. We then ex- 
amined properties of genotypes every n steps of this ran- 
dom walk, where n equaled 5 times the upper limit m + 
of the number of interactions in the circuit. This sporadic 
sampling serves to erase correlations in genotypes along 
this random walk. 
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Results 



dynamics, and lead to a new steady-state activity pattern 



Genotype networks of gene expression phenotypes 
have different sizes 

For our model, most or all genotypes that produce the 
same p henotype form large connected genotype net- 
works (Cilibert i et all l2Q07alTbb . The size of any one 
phenotype's genotype network depends only on the frac- 
tion d of genes whose expression state differs between 
the initial s tate Sp, and the stead y state activity pheno- 
type s££* dCiliberti et all l2007bl) . Specifically, pheno- 
types where these two states (regardless of their actual 
expression values) are more similar have larger genotype 
networks (Fig. SI in Supplementary Material). One can 
view regulatory circuits as devices that compute an ex- 
pression state s°§* from the initial state so. From this 
perspective, a larger number of gene expression differ- 
ences between these states means that the computation 
becomes increasingly difficult, in the sense that fewer 
genotypes can perform it. 

We first examined, in our model, the relationship 
between the size of a phenotype P's genotype network 
and the robustness of circuits with this phenotype P to 
mutations. To this end, we uniformly sampled 10 6 geno- 
types from genotype networks of different sizes (differ- 
ent d), and determined their mean robustness to muta- 
tions, that is, the mean fraction of their neighbors with 
the same phenotype. For all examined cases, the aver- 
age mutational robustness (i.e. phenotypic robustness) is 
higher for genotypes on larger genotype networks (Fig. 
S2). Thus, phenotypic robustness to mutations in creases 
with g enotype network size, just as for RNA dWagneil 
2008b). Therefore, we can simply use 1 — d as a proxy 



for genotype network size and phenotypic robustness to 
mutations. 

Phenotypic robustness to mutations facilitates pheno- 
typic variability in response to noise 

In this paper, we are concerned with the production 
of new steady-state gene expression patterns s^ w that 
are different from s££*. We refer to such activity pat- 
terns as new phenotypes. They could result from mu- 
tations that change regulatory interactions in a circuit. 
They could also result from non-genetic perturbations 
(Fig. lc). We here consider two kinds of non-genetic per- 
turbations, noise in a cell's internal milieu, and change in 
the organism's (external) environment. Both kinds can 
induce dramatic gene expression c hanges in organisms 
ranging from bacteria to metazoans (jElowitz et all l2002t 



Rai etall l20ld: ISnell-Rood etail [2010b. We first focus 



on noise, which includes stochastic changes in protein or 
mRNA copy numbers in a cell, and can cause phenotypic 
heterogeneity in clonal populations jFTowitz et al , 2002 ; 



iMcAdams "and ArkinLll997tlRaj et alll2010h . Such noise 
may affect the activity or expression of circuit genes at a 
given time, which may alter a circuit's gene expression 



We emulated the perturbations produced by noise in 
two complementary ways. First, we changed the activity 
state of single genes in the initial state sq> f° r eacn gene 
in a circuit, and determined the different new phenotypes 
g new t j lat resu it e( j f rom suc h change. Secondly, we ran- 
domly perturbed the dynamic trajectory from so to Sqq 
('noisy dynamics'), as described in Methods. 

We asked how the mutational robustness of a gene 
expression phenotype affects the number of new pheno- 
types that these two kinds of noise can produce in pop- 
ulations of evolving circuits. We evolved populations of 
200 circuits under stabilizing selection on a given gene 
expression state s^*, as described in Methods. We then 
counted the number of unique new phenotypes that the 
two different kinds of noise produced among the individ- 
uals in a population. Specifically, we counted phenotypes 
that appeared multiple times only once, thus focusing on 
unique phenotypes. 

Noise can produce a greater number of new phe- 
notypes in populations evolving on large genotype net- 
works. Fig. 2 demonstrates these observations for cir- 
cuits with N = 20 genes and an interaction density 
c « 0.2. These observations also hold if we vary the 
numbers of genes and regulatory interactions in a circuit 
(Figs. S3,4), with a single exception for perturbations 
in so when the number of regulatory interactions is very 
low (Fig. S3d). 

Populations with more robust phenotypes harbor 
more diverse genotypes 

Increased genotypic diversity in populations evolving 
in large genotype networks might aid in producing in- 
creased phenotypic variability, as discussed in the Intro- 
duction. We next asked whether this mechanism may 
apply to our system. To this end, we quantified the geno- 
typic distance among two circuits (see Methods). 

As a measure of a population's genotypic diversity, 
we estimated the mean pairwise circuit distance, as well 
as its maximum, in each of 500 populations evolved un- 
der stabilizing selection on a phenotype s££*. We did so 
for two classes of populations that differ in the robustness 
of their phenotypes, and found that the mean genotypic 
distance is significantly higher for populations with a ro- 
bust phenotype. The same holds also for the maximum 
genotypic distance. These observations are not sensitive 
to the number of genes and interactions in a circuit (Ta- 
ble SI). Thus, populations with a robust phenotype are 
genetically more diverse than populations with a less ro- 
bust phenotype. These observations hint that the higher 
genetic diversity of populations with robust phenotypes 
may be exposed as phenotypic variability in response to 
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Fig. 2. High phenotypic robustness facilitates phenotypic variability in 
response to noise in gene expression. The distance d between the 
initial state so and the optimal phenotype s£o* is strongly associated 
with genotype network size and with phenotypic mutational 
robustness ("phenotypic robustness") hereafter. 'High', 'medium' and 
'low' correspond to expression phenotypes with high (d = 0.1), 
intermediate (d = 0.25), and low (d = 0.5) robustness. The figure 
shows results for N = 20 genes and a fraction c ss 0.2 of non-zero 
regulatory interactions. Both panels show mean numbers of novel 
phenotypes averaged over 500 independent populations, for each level 
of robustness. The length of bars denotes one standard error. The 
number of different new phenotypes that a population can access after 
perturbations of a) single genes in the initial state so, or b) a circuit's 
gene expression trajectory, increases with phenotypic robustness. 



Phenotypic robustness does not facilitate phenotypic 
variability caused by mutations 



We next asked whether phenotypic robustness also facil- 
itates phenotypic variability in response to mutations for 
the regulatory circuits we study. We again studied popu- 
lations of circuits evolved under stabilizing selection on 
a phenotype s°£*. In such populations, we determined 
the 1 -mutation neighborhood of each circuit, that is, all 
circuits that differ from it in a single regulatory interac- 
tion. We then determined the number of unique new gene 
expression phenotypes in the population's neighborhood. 
That is, we counted only once a phenotype if it occurred 
in the neighborhood of two different circuits. This num- 
ber of unique phenotypes is a measure of the population's 
phenotypic variability in response to mutations. 
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Fig. 3. High phenotypic robustness does not facilitate phenotypic 
variability in response to mutations without preceding environmental 
change. 'High', 'medium' and 'low' correspond to expression 
phenotypes with high (d = 0.1), intermediate (d = 0.25), and low 
(d = 0.5) robustness. The figure shows results for N = 20 genes and 
a fraction c m 0.2 of non-zero regulatory interactions. The panel 
shows the mean number of novel phenotypes averaged over 500 
independent populations, at each level of robustness. The length of 
bars denotes one standard error. 

We found that populations with a highly robust phe- 
notype show lower phenotypic variability in response to 
mutations. This holds despite their somewhat higher 
genotypic diversity (Table SI, discussed above). Fig. 3 
shows pertinent data for circuits with 20 genes and inter- 
action density c « 0.2. The same behavior holds for 
populations of circuits with different number of genes 
and different interaction densities (Fig. S5). In sum, ro- 
bustness of a phenotype to mutations impairs phenotypic 
variability to mutation, as opposed to what we saw for 
variability in response to noise. 

Our results suggest that a phenotype's mutational 
robustness promotes phenotypic variability in response 
to noise, but hinders such variability in response to mu- 
tations. This may seem surprising, because robustness 
to mutations inc reases with robustnes s to noise for indi- 
vidual circuits (Ciliberti et al. . l2007bl) . One might thus 
think that phenotypic variability also behaves similarly 
in response to these perturbations. However, robustness 
to mutations explains less than 25 percent of the vari- 
ance in robustness to noise, as shown by a new statistical 
analys es of our previously published data dCiliberti et all 
2007bJ). Thus phenotypic variability in response to noise 
and to mutation are only weakly coupled. 

With these observations in mind, we analyzed the 
phenotypic variability in response to noise and mutations 
of individual circuits in populations evolving on differ- 
ent genotype networks (Table S2). After obtaining this 
data, we compared the mean number of new phenotypes 
that mutations or noise could produce from circuits in 
populations with different levels of phenotypic robust- 
ness (Table S3). We found that phenotypic variability 
in response to gene expression noise decreases less with 
phenotypic robustness to mutations than phenotypic vari- 
ability to mutations (Table S3). It may even increase with 



phenotypic robustness. These observations suggest that 
the increased genotypic diversity attained on larger geno- 
type networks is insufficient to compensate for the small 
(or null) reduction in variability in response to mutations. 
It is, however, sufficient to compensate for the smaller re- 
duction in phenotypic variability in response to noise in 
gene expression. 

Phenotypic robustness increases phenotypic variabil- 
ity after environmental change 

Thus far, we focused mostly on phenotypic variability 
in response to small, random non-genetic perturbations, 
such as single gene expression perturbations along a gene 
expression trajectory. We now turn to the question of 
what happens when a whole population is subject to the 
same non-genetic perturbation. In nature, this may occur 
because of environmental change outside the organism, 
or colonization of a new habitat. 

The environment can have two different roles in this 
context. The first is an inducing role, w here the environ- 
ment acts as an "agent of development' ' dWest -Eberhard, 
1989). In this role, it affects the phenotype produced 
from a genotype. In many cases, environmentally in- 
duced phenotypi c change is linked to ma jor changes in 
gene expression (iSnell-Rood et all 120101) . The second 
role is an evaluating role, where the envi r onme nt acts as 
an "agent of selection" ( West-Eberhard , 1989h . In an- 
thropomorphic terms, the environment in this role distin- 
guishes well-adapted from poorly adapted phenotypes. 

Conveniently, our model allows us to study these 
roles independently. We model a change in the environ- 
ment's evaluation role as a change in the identity of the 
optimal phenotype s^*, for all circuits in the population. 
We model a change in the environment's inducing role 
as a change in the initial state sq in the whole popula- 
tion. Such a change could occur, for example, through a 
signaling pathway that detects an environmental change, 
and that affects genes upstream of the circuit. Put dif- 
ferently, changes in sq reflect the environment's effect 
on phenotype production, while changes in s°£* affect 
the survival probability of individuals, without inducing 
novel phenotypes. We note that other factors, such as 
mutations in upstream genes, might also lead to changes 
in so- Any one such change, however, would initially 
affect only one individual in a population, and not the 
whole population at the same time. 

We first asked how an environmentally induced 
change in the initial gene activity pattern so affects the 
number of different actual phenotypes that a population 
displays. We note that our populations may contain a 
few individuals with phenotypes different from the op- 
timal phenotype s^*. The reason is that, in contrast to 



der stabilizing selection on an optimal expression pheno- 
type s°g* and a given gene activity pattern sg as initial 
condition. We counted the number of phenotypes in the 
population, and compared it with the number of different 
phenotypes that the same population displays when sg is 
replaced by a random gene activity pattern sg as an initial 
condition. We found that phenotypic variability increases 
after substitution of Sq with Sq (Figs. 4 and S6). In addi- 
tion, the magnitude of this increment increases with phe- 
notypic robustness (Fig. 4). This last observation is gen- 
erally not sensitive to the number of genes and regulatory 
interactions in a circuit (Fig. S6). The increase in phe- 
notypic variability is significantly higher for populations 
with a robust phenotype (Mann- Whitney U-test; Table 
S4). The single exception to these observations were cir- 
cuits of very low interaction density (N = 20; c ss 0.1 
Fig. S6d), tha t also s how other non-typical behaviors 



(ICiliberti et all 12007a) . Our results suggest that, after 



previous formulations dCiliberti et al. , 2007b ). we here 
represent fitness as a continuous variable that depends 
on the similarity of a circuit's phenotype Soo to sJ^J* (see 
Methods). We started out with a population evolved un- 



environmental change, observable phenotypic diversity 
increases to a larger extent in populations with a robust 
phenotype. We note that because the identity of does 
not affect the production of phenotypes, but only their 
viability, it is not appropriate to carry out an analogous 
analysis for changes in s^J*. 

In earlier sections, we have shown that phenotypic 
robustness impedes phenotypic variability after muta- 
tions in populations evolving in a constant environment 
(Fig. 3). We next asked whether this also holds after 
a change in the inducing role of the environment. We 
started out, as in our last analysis, with populations of 
circuits evolved under stabilizing selection on an optimal 
expression phenotype and with a given gene activ- 
ity pattern sg as initial condition. Then, we changed the 
initial condition sg for all the circuits to a new random 
initial condition Sq, and allowed evolution to proceed. 
Before and after this change, we recorded the number 
of different phenotypes accessible from the population 
through mutations. Under the new condition the popula- 
tion effectively searches genotype space for optimal phe- 
notypes. During this search, many variant circuits may 
not pass to subsequent generations. Our primary focus, 
however, is not this search, but the immediate increase in 
a population's phenotypic variability in response to envi- 
ronmental change. 

Before environmental change, populations with a 
robust phenotype have access to fewer phenotypic vari- 
ants, just as in our previous observations (Fig. 3). Imme- 
diately after environmental change (at t = 1), the num- 
ber of new phenotypes accessible through mutations in- 
creases, in a burst, in all populations. Importantly, this 
increase is higher in populations with a robust phenotype 
(Figs. 5a and S7). This means that phenotypic robust- 
ness facilitates the phenotypic variability caused by mu- 
tations, but only after environmental change. As in our 
analysis above, the only exception occurs when interac- 
tion density is very low (Fig. S7d). 
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low medium high 

Phenotypic robustness 

Fig. 4. High phenotypic robustness increases phenotypic diversity in 
populations of gene circuits after environmental change. The number 
of different phenotypes that populations display increases after 
changing the initial expression state so- Such an increase is greater for 
populations with mutationally more robust phenotypes. The figure 
shows results for N = 20 genes and a fraction c 0.2 of non-zero 
regulatory interactions. The panel shows the mean number of observed 
phenotypes averaged over 500 independent populations, for each level 
of robustness. The length of bars denotes one standard error. 

We next asked whether the evaluation role of the en- 
vironment has similar effects on the mutational access to 
new phenotypes. To this end, we repeated the above anal- 
ysis, but replaced, at t = 1, the optimal phenotype s^ t,a 
by a randomly chosen optimal s°^ t,b (without changing 
so). We also observed a transient, but more gradual, in- 
crease in the number of phenotypes that are mutationally 
accessible. However, in this case, phenotypic variabil- 
ity is lower for populations with a robust phenotype after 
environmental change (Figs. 5b and S8). Thus, the in- 
ductive role of environment, but not its evaluative role, 
causes higher phenotypic variability in populations with 
robust phenotypes. 

An open question is how mutation-accessible phe- 
notypic variability changes when both the inductive and 
the evaluation roles of the environment change. This 
question is important because a change in the inductive 
role favors phenotypic variability to a larger extent in 
populations with a robust phenotype, whereas a change 
in the evaluative role particularly favors variability in 
populations with less robust phenotypes. Thus, a combi- 
nation of both effects could result in a negligible effect of 
phenotypic robustness on variability after environmental 
change. To answer this question, we repeated our analy- 
sis from the previous paragraph, but replaced the original 
pair of states (sg, s°^ t,a ) with a new pair (s\, s^*' b ), such 
that the distance d between sq and s^* was the same for 
both pairs. 

In this new analysis, populations evolving on a large 
genotype network show greater phenotypic variability in 
response to mutations immediately after this change (Fig. 
5c). These differences are statistically highly significant 
(Table S5). The same observations hold for circuits of 
different sizes and different connectivities (Fig. S9 and 
Table S5). As in our analysis above, the only exception 



occurs when interaction density is very low (Fig. S9d). 
These observations imply that the inductive role domi- 
nates in its immediate effect on phenotypic variability 
when both roles of the environment change. In sum, pop- 
ulations with a robust phenotype have mutational access 
to more phenotypic variants after environmental change. 
This increased access is caused by the inductive role of 
the environment, that is by the new phenotypes that a new 
environment can bring forth. 

Discussion 

If non-genetic change is to be causally involved in evo- 
lutionary innovation, it needs to generate novel phe- 
notypes. Genetic assimilation can then stabilize these 
phenotypes if the non-genetic perturbations that induced 
them appear recurrently. Here, we focused on whether 
the robustness of an existing phenotype and non-genetic 
change can facilitate the origin of new phenotypes. To 
address this question, we examined a generic model of 
transcriptional regulation circuitry, in which the rela- 
tionship between genotypes (patterns of regulatory in- 



sion patterns) is well-studied ( 


Ciliberti et all 


2007a Ibt 


Martin and Wasner, 2008, 2009 


; Waenei. 1996 


). In this 



model, we can use the size of a phenotype's genotype 
network as a proxy for a phenotype's robustness to mu- 
tations. 

We broadly distinguished two kinds of non-genetic 
perturbations. The first corresponds to fluctuations in 
a gene circuit's microenvironment that produce random 
changes in gene expression. Such changes, are ubiqui- 
tous an d have important effects on cell biological pro- 
cesses (Elowitz et ail 120021: iMcAdams and ArkinLll997t 



Raj et all |2010). The second kind comprises changes 
in the (macro)environment external to an organism. For 
brevity, we refer to these kinds of change as noise and 
environmental change. 

We first explored how noise and mutations affect 
phenotypic variability in populations that differ in the 
robustness of their gene activity phenotypes. We found 
that phenotypic mutational robustness increases pheno- 
typic variability of populations in response to noise but 
not in response to mutations. This last finding differs 
from observations for RNA secondary structure, where 
phenotypic robustness f acilitates the m utational access 
to phenotypic variants (IWagnerl l2008bl) . The reason 
stems from differences in the organization of genotype 
space for these two system classes, i.e. in the dis- 
tributi on of different genotype networks in genotype 
space dAncel and Fontanal 2000; Ciliber ti et all , 2007a; 



Espinosa-Soto et all 1 20 10t IWagnerl I2008bl) 7 A recent 



mathematical model (Dra ghi et all 12010 ) shows that mu- 
tational access to new phenotypes can depend on this or- 
ganization, and on details of a population's evolutionary 
dynamics. 

Next, we asked how phenotypic robustness af- 
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Fig. 5. High phenotypic robustness allows mutational access to more phenotypes after an environmental change produces novel phenotypes.The 
figure shows the number of different phenotypes in the 1 -mutant neighborhood of a population, for three different scenarios of environmental 
change at generation t = 1. The insets show the number of phenotypes accessible through mutation immediately before (t = 0) and 
immediately after (t = 1) environmental change. The figure shows results for N = 20 genes and a fraction c 0.2 of non-zero regulatory 
interactions. 'High', 'medium' and 'low' correspond to expression phenotypes with high (d = 0.1), intermediate (d = 0.25), and low (d = 0.5) 
phenotypic robustness. Data are mean values averaged across 500 independent simulations for each level of robustness. The length of bars 
denotes one standard error, (a) The initial state so changes, (b) The identity of the optimal phenotype s^J* changes, (c) The original pair (Sq, 



s oo ,a ) changes to a new pair (sjj, 



opt,b 



fects phenotypic variability in response to environmen- 
tal change, which we modeled to affect all individuals 
in a population. We showed that environmental change, 
besides increasing observable phenotypic variation, tran- 
siently increases phenotypic variability caused by mu- 
tations. Because mutational access to most novel vari- 
ants is only possible in the new environment, these vari- 
ants can be considered environmentally-induced pheno- 
types, supporting the phenotype-first scenario. From this 
perspective, environmental change increases phenotypic 
variability and the chances to refine, exaggerate and sta- 
bilize phenotypic variation via genetic change. Impor- 
tantly, this increase in phenotypic variability is higher in 
populations that had a more robust phenotype before en- 
vironmental change. 

In sum, we found a positive effect of phenotypic 
robustness on phenotypic variability after non-genetic 
perturbations. This positive effect is not sensitive to 
the magnitude of non-genetic perturbations. Phenotypic 
robustness favors phenotypic variability after single- 
gene perturbations in the initial condition, which is the 
smallest possible non-genetic perturbation in our model 
(Fig. 2a). It also favors phenotypic variability after re- 
placing the initial condition by a completely new ran- 
dom activity pattern (Figs. 4 and 5a). In contrast, pop- 
ulations with robust gene expression phenotypes pro- 
duce less phenotypic variation when only mutations are 
used to explore the phenotypic possibilities. Thus, a 
mechanism that relies exclusively on mutation to pro- 
duce novel phenotypes becomes less important for inno- 
vation as a phenotype's robustness increases. Our results 
suggest that plasticity-mediated innovation may be espe- 
cially important for gene expression traits with high mu- 
tational robustness. Our work is an initial step towards 
the definition of two different domains where either the 
genotype-first or phenotype-first scenarios prevail over 
the other. In this regard, we note that environmental in- 
duction of novel traits is not only expected for gene cir- 



cuits with robust phenotypes. In fact, we observe a gen- 
eral increase in phenotypic variability after environmen- 
tal change (Figs. 4 and 5). It is just that this increase is 
especially strong for robust phenotypes. 

The general increase in phenotypic variability 
we observe is consistent with many empirical ob- 
servations on phenotypic variation that is conditional 
on the environment. For example, severe environ- 
ments e nhance phenotypic diff e rence s among fruit fly 
strains (IKondrashov and Houlel 1 19941) . and a temper- 
ature rise due to an unshaded milieu increases the 
frequency of abnormal morphologies in fruit flies 
(Roberts and Feder, Il999h . Vloreover, populat on genet- 
ics studies show that the release of hidden genetic varia- 
tion after environmental chang e should be very common 
dHermisson and Wagneru2004l) . 

Unfortunately, because the mutational robustness 
of most traits is unk nown and difficult to measure 
(Ide Visser et al ., 2003), these observations do not speak 
directly to the question whether phenotypic robust- 
ness promotes phenotypic variability after environmen- 
tal change. However, there is an intriguing phenomenon 
with potential parallels to our observations. In the 
fruit fly D. melanogaster and in the plant Arabidopsis 
thaliana, previously unobserved phenotypic variants ap- 
pear when the activity of the chaperone protein Hsp90 
is impaired (by environmental stress or by other means ) 



dOueitsch et all 2002; Ruther ford and Lindquist 119981) . 



This observation suggests that an active Hsp90 protein 
usually conceals genetic variation that would otherwise 
be visible phenotypically, although the appearance of 
new genetic varia nts may underlie some of the new phe- 
notypic variation dSpecchia et al. ■ l2010l) . Hence, Hsp90 
might increase phenotypic variability by allowing the ac- 
cumulation of cryptic genetic variation that can become 
phenotypically visible in the right environment. Our 
analysis revealed an analogous phenomenon for robust 
gene expression phenotypes. Robustness of phenotypic 
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traits can facilitate the accumulation of cryptic genetic 
variation in a population. Non-genetic perturbations can 
expose this variation as new phenotypes, some of which 
may be adaptive. The parallels between the action of 
Hsp90 and the phenomenon we study may be more than 
superficial. The reason is that Hsp90 increases the num- 
ber of genotypes w i th a given phenotype for several t raits 
feufherfordl. 120031 iRutherford and Lindquistl Il998l) . It 
thus increases phenotypic robustness. Hsp90 is a specific 
example - with a peculiar mechanism - of how pheno- 
typic robustness can increase phenotypic variability after 
non-genetic perturbations. 

In conclusion, our observations suggest that pheno- 
typic robustness to mutations can play a positive role 
in the phenotypic variability after non-genetic pertur- 
bations. To see this, one needs to study the role of 
population level processes, as we did. We caution that 
we made our observation in the context of a specific 
model of transcriptional regulation circuits. The gene 
expression phenotypes of such circ uits play central roles 
in many evolutionary inn ovations ([Davidson and Erwin , 
2006; IShubin et al.[|2009l) . However, phenotypes may be 
distributed differently in genotype space in other classes 
of biological systems. Thus, whether our observations 
hold in other systems remains to be seen. 
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